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Abstract. Mass-action chemical reaction systems are frequently used in Computational Biology. 
The corresponding polynomial dynamical systems are often large (consisting of tens or even hun- 
dreds of ordinary differential equations) and poorly parametrized (due to noisy measurement data 
and a small number of data points and repetitions). Therefore, it is often difficult to establish the 
existence of (positive) steady states or to determine whether more complicated phenomena such as 
multistationarity exist. If, however, the steady state ideal of the system is a binomial ideal, then 
we show that these questions can be answered easily. The focus of this work is on systems with this 
property, and we say that such systems have toric steady states. Our main result gives sufficient 
conditions for a chemical reaction system to have toric steady states. Furthermore, we analyze the 
capacity of such a system to exhibit positive steady states and multistationarity. Examples of sys- 
tems with toric steady states include weakly-reversible zero-deficiency chemical reaction systems. 
An important application of our work concerns the networks that describe the multisite phospho- 
rylation of a protein by a kinase/phosphatase pair in a sequential and distributive mechanism. 

Keywords: chemical reaction networks, mass-action kinetics, multistationarity, multisite phospho- 
rylation, binomial ideal. 
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1. Introduction 



Ordinary differential equations (ODEs) are an important modeling tool in Systems Biology and 
many other areas of Computational Biology. Due to the inherent complexity of biological sys- 
tems, realistic models are often large, both in terms of the number of states and the (unknown) 
parameters. Moreover, models are often poorly parametrized, a consequence of noisy measurement 
data, a small number of data points, and a limited number of repetitions. Hence, for mass-action 
chemical reaction systems, the focus of the present article, simply establishing the existence of 
(positive) steady states can be demanding, as it requires the solution of a large polynomial system 
with unknown coefficients (usually the parameters). Moreover, due to the predominant parameter 
uncertainty, one is often not interested in establishing the existence of a particular steady state, 
but rather in obtaining a parametrization of all steady states - preferably in terms of the sys- 
tem parameters [34J. Frequently one is also interested in the existence of multiple steady states 
(multistationarity), for example, in modeling the cell cycle [21 [31 [28], signal transduction [23[ [26] 
or cellular differentiation |32[ I33j. For general polynomial systems with unknown coefficients, the 
tasks of obtaining positive solutions or a parametrization of positive solutions, and deciding about 
multiple positive solutions, are clearly challenging. For the systems considered in this article - 
chemical reaction systems with toric steady states - these questions can be answered easily. 

We say that a polynomial dynamical system dx/dt = f(x) has toric steady states if the ideal 
generated by its steady state equations is a binomial ideal (see Definition 12. 2p . We give sufficient 
conditions for a chemical reaction system to have toric steady states (Theorems 13.81 and I3.19P and 
show in this case that the steady state locus has a nice monomial parametrization (Theorems 13.111 
and 13.20]) . Furthermore, we show that the existence of positive steady states in this case is straight- 
forward to check (Theorem I5.5p . 

There are several important classes of mass-action kinetics chemical reaction systems which have 
toric steady states. These include usual instances of detailed-balanced systems in the sense of 
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Feinberg, Horn, and Jackson fT2j [20l I21j . which show particularly nice dynamical behavior. 
These systems are weakly-reversible, a hypothesis we do not impose here. 

A chemical reaction system with toric steady states of great biological importance is the mul- 
tisite phosphorylation system; this network describes the n-site phosphorylation of a protein by a 
kinase/phosphatase pair in a sequential and distributive mechanism. Biochemically, these systems 
play an important role in signal transduction networks, cell cycle control, or cellular differentiation: 
for example, members of the family of mitogen-activated kinase cascades consist of several such 
phosphorylation systems with n = 2 or n = 3 (see e.g. [221129]). the progression from Gl to S phase 
in the cell cycle of budding yeast is controlled by a system with n = 9 (by way of the protein Sicl, 
see e.g. [9]), and a system with n = 13 plays an important role in T-cell differentiation (by way of 
the protein NFAT [171 HSJ |22|). 

Consequently there exists a body of work on the mathematics of phosphorylation systems and 
the more general class of post-translational modification systems: for example, Conradi et al. [6], 
Wang and Sontag [36], Manrai and Gunawardena [25], and Thomson and Gunawardena |34[ 135]. 
While the first two references are concerned with the number of steady states and multistationar- 
ity, the references of Gunawardena et al. deal with parametrizing all positive steady states. The 
present article builds on these earlier results. In fact, the family of monomial parametrizations 
obtained here for multisite phosphorylation systems (Theorem 14. 3|) is a specific instance of a ra- 
tional parametrization theorem due to Thomson and Gunawardena, and one parametrization of 
the family was analyzed earlier by Wang and Sontag. Furthermore, we show that by using results 
from [6] one can determine whether multistationarity exists for systems with toric steady states by 
analyzing certain linear inequality systems. In this sense our results can be seen as a generalization 
of @. 

This article is organized as follows. Section [2] provides an introduction to the mathematics of 
chemical reaction systems. Our main results on toric steady states appear in Section [3l Theo- 
rems 13.81 and 13.191 give sufficient criteria for a system to exhibit toric steady states, and Theo- 
rems [3JJ] and [3]20] give parametrizations for the steady state locus. As an application of this work, 
we analyze the steady state loci of multisite phosphorylation systems in Section HI Theorem 14.31 
summarizes our results: we show that these systems have toric steady states for any choice of re- 
action rate constants, and we give an explicit parametrization of the steady state locus. Section [5] 
focuses on multiple steady states for chemical reaction systems with toric steady states. Theo- 
rem [53] gives a criterion for such a system to exhibit multistationarity, and we make the connection 
to a related criterion due to Feinberg. 

2. Chemical reaction network theory 

In this section we recall the basic setup of chemical reaction systems, and we introduce in § 12.21 
the precise definition of systems with toric steady states. We first present an intuitive example 
that illustrates how a chemical reaction network gives rise to a dynamical system. An example of 
a chemical reaction, as it usually appears in the literature, is the following: 

A + B ^^3A + C (2.1) 

In this reaction, one unit of chemical species A and one of B react (at reaction rate k) to form 
three units of A and one of C. The educt (or reactant or source) A + B and the product 3A + C are 
called complexes. We will refer to complexes such as A + B that are the educt of a reaction as educt 
complexes. The concentrations of the three species, denoted by xa, xb, and xc, will change in time 
as the reaction occurs. Under the assumption of mass- action kinetics, species A and B react at a 
rate proportional to the product of their concentrations, where the proportionality constant is the 
rate constant k. Noting that the reaction yields a net change of two units in the amount of A, we 
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obtain the first differential equation in the following system: 



d 
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The other two equations arise similarly. A chemical reaction network consists of finitely many 
reactions. The differential equations that a network defines are comprised of a sum of the monomial 
contribution from the reactant of each chemical reaction in the network; these differential equations 
will be defined in equation (12. 3p . 

2.1. Chemical reaction systems. We now provide precise definitions. A chemical reaction net- 
work is a finite directed graph whose vertices are labeled by complexes and whose edges are labeled 
by parameters (reaction rate constants). Specifically, the digraph is denoted G = (V,E), with 
vertex set V = {1, 2, ... , m} and edge set E C {(i,j) £ V x V : i ^ j}. Throughout this article, 
the integer unknowns m, s, and r denote the numbers of complexes, species, and edges (reactions), 
respectively. Linkage classes refer to the connected components of a network, and terminal strong 
linkage classes refer to the maximal strongly connected subgraphs in which there are no edges 
(reactions) from a complex in the subgraph to a complex outside the subgraph. The vertex i of G 
represents the i-th chemical complex, and we associate to it the monomial 



More precisely, if the i-th complex is y%\A + y^E + • • • (where E Z>o for j = 1, 2, . . . , s), then 
it defines the monomial x A ll x y ^ 2 ■ ■ ■ . For example, the two complexes in the network (12. ip give rise 
to the monomials xaxb and x\xc, which determine two vectors y\ = (1, 1,0) and 7/2 = (3,0, 1). 
These vectors define the rows of an m x s-matrix of non- negative integers, which we denote by 
y = (Uij)- Next, the unknowns x\,X2, ■ ■ ■ ,x s represent the concentrations of the s species in the 
network, and we regard them as functions Xi(t) of time t. The monomial labels form the entries in 
the following vector: 



A directed edge £ E represents a reaction from the i-th chemical complex to the j'-th 

chemical complex. Each edge is labeled by a positive parameter Kij which represents the rate 
constant of the reaction. In this article, we will treat the rate constants as unknowns; we are 
interested in the family of dynamical systems that arise from a given network as the rate constants 
vary. 

The main application of our results are chemical reaction networks under mass- action kinetics. 
Therefore, even if the principal results in § [3] hold for general polynomial dynamical systems, we 
assume in what follows mass-action kinetics. We now explain how mass-action kinetics defines a 
dynamical system from a chemical reaction network. Let A K denote the negative of the Laplacian 
of the chemical reaction network G. In other words, A K is the m x m-matrix whose off-diagonal 
entries are the Kij and whose row sums are zero. Now we define the complex-to-species rate matrix 
of size s x m to be 



x 




■ Vis 




S := 





The reaction network G defines the following dynamical system: 

dx fdxidx2 dx s V 

dt \ dt ' dt dt J ' ^'^^l 



(2.3) 
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We see that the right-hand side of each differential equation dxi/dt is a polynomial in the polynomial 
ring M[(Ky)(jj) g £, xi,x 2 , . . . ,x s }. A chemical reaction system refers to the dynamical system (I2.3P 
arising from a specific chemical reaction network G and a choice of rate parameters (n*j) € M> 
(recall that r denotes the number of reactions). 



Example 2.1. The following chemical reaction network is the 1-site phosphorylation system: 



S + E 

k 



onn 



ESo kc 4" Si + E 



(2.4) 



offo 



S + F . 



The key players in this network are a kinase enzyme (E), a phosphatase enzyme (F), and two 
substrates (So and Si). The substrate S\ is obtained from the unphosphorylated protein So by 
attaching a phosphate group to it via an enzymatic reaction involving E. Conversely, a reac- 
tion involving F removes the phosphate group from S\ to obtain So. The intermediate com- 
plexes ESo an d ES± are the bound enzyme- substrate complexes. Under the ordering of the 6 species 
as (So, S\, ESo, FSi, E, F) and the 6 complexes as (So + E, Si + E, ESo, So + F, Si + F, FSi), the 
matrices whose product defines the dynamical system \2. 3\) follow: 



^(x) 



(x So X E , X Sl X E , XES , x So x f , x Sl x F , X FSl ] 
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We will study generalizations of this network in this article. 

The stoichiometric subspace is the vector subspace spanned by the reaction vectors yj —yi (where 
(i,j) is an edge of G), and we will denote this space by S: 

S := M{ yj - yi | (i,j)GE} . 

In the earlier example shown in (|2.ip . we have y 2 — yi = (2, —1, 1), which means that with the 
occurrence of each reaction, two units of A and one of C are produced, while one unit of B is 
consumed. This vector (2,-1,1) spans the stoichiometric subspace S for the network (12. ip . Note 



that the vector ^ in (|2.3p lies in S for all time t. In fact, a trajectory x(t) beginning at a positive 
vector x(0) = x° G W >0 remains in the stoichiometric compatibility class (also called an "invariant 
polyhedron"), which we denote by 



(x° + s) n : 



s 

^>0 > 



(2.5) 



for all positive time. In other words, this set is forward-invariant with respect to the dynamics f|2.3[) . 
It follows that any stoichiometric compatibility class of a network has the same dimension as the 
stoichiometric subspace. 
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2.2. Steady states. We present the definition of systems with toric steady states. For background 
information on the algebraic tools we use, we refer the reader to the nice textbook of Cox, Little, 
and O'Shea 0. 

Recall that an ideal in M[a;i, x 2 , ■ ■ ■ , x s ] is called a binomial ideal if it can be generated by 
binomials (i.e., polynomials with at most two terms). The basic building blocks of binomial ideals 
are the prime binomial ideals, which are called toric ideals |10| . 

Definition 2.2. Consider a polynomial dynamical system dxi/dt = fi(x), for i = 1, 2, . . . , s, with 
fli /2> • • • j fs £ M.[x\,X2, ■ ■ ■ , x s ]. We are interested in the real zeros of the steady state ideal: 

Jy,v = (fi,f2, ■■■,/ s )= 71 hj{x)fj(x) | hi(x) € R[xi,x 2 , ... ,x s ] for 1 < % < s 

U=i 

The real zeros of Js* are called steady states, and the term steady state locus is used to denote 
the set of real zeros of J^q, : 

{x*GR s | fi(x*) = f 2 (x*) = --- = f e (x*) = 0}. 

We say that the polynomial dynamical system has toric steady states if Js<i> is a binomial ideal and 
it admits real zeros. 

We are interested in positive steady states x £ R^q and will not be concerned with boundary 
steady states x G (IR> \ M> ) . 

This article focuses on mass-action kinetics chemical reaction systems. In this case, the poly- 
nomials fi, f 2 , ■ ■ ■ , f s correspond to the rows of the system (12.31) . In general, having toric steady 
states depends both on the reaction network and on the particular rate constants, as the following 
simple example shows. 

Example 2.3 (Triangle network). Let s = 2, m = 3, and let G be the following network: 



2A 




K32 



We label the three complexes as x yi = x\, x y2 = x\, x vz = x\x 2 , and we define Kij to be the (real 
positive) rate constant of the reaction from complex x Vi to complex x y K The resulting mass-action 
kinetics system (|2.3p equals 

~d~t = ~ = (~ 2ki2 ~ Kl 3) x i + ( 2k 21 + ^23)^ + («3i - k 32 )xix 2 • 

Then, the steady state locus in M? is defined by this single trinomial. As only the coefficient of x\x 2 
can be zero, this system has toric steady states if and only if K31 = K32 . 

A chemical reaction system exhibits multistationarity if there exists a stoichiometric compat- 
ibility class V x o with two or more steady states in its relative interior. A system may admit 
multistationarity for all, some, or no choices of positive rate constants Kij] if such rate constants 
exist, then we say that the network has the capacity for multistationarity. 

2.3. The deficiency of a chemical reaction network. The deficiency S of a chemical reaction 
network is an important invariant. For a chemical reaction network, recall that m denotes the 
number of complexes. Denote by I the number of linkage classes. Most of the networks considered 
in this article have the property that each linkage class contains a unique terminal strong linkage 
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class. In this case, Feinberg showed that the deficiency of the network can be computed in the 
following way: 

5 := m — I — dim(iS) , 

where S denotes the stoichiometric subspace. Note that in this case the deficiency depends only 
on the reaction network and not on the specific values of the rate constants. The deficiency of a 
reaction network is non-negative because it can be interpreted as the dimension of a certain linear 
subspace or the codimension of a certain ideal [8]. For systems arising from zero-deficiency 
networks and networks whose linkage classes have deficiencies zero or one, there are many results 
due to Feinberg that concern the existence, uniqueness, and stability of steady states [TT | \12 \ [TB I fT3] . 

3. Sufficient conditions for the existence of toric steady states 

The main results of this section, Theorems 13.31 13. 8\ and I3.19| give sufficient conditions for a 
chemical reaction system to have toric steady states and state criteria for these systems to have 
positive toric steady states. Theorems 13.111 and 13 . 201 give a monomial parametrization of the steady 
state locus in this case. 

We first state several conditions and intermediate results that will lead to Theorem 13.81 Recall 
that a partition of {l,2,...,m} is a collection of nonempty disjoint subsets I\, I2, ■ ■ ■ , Id with 
respective cardinalities l\, I2, ■ ■ ■ , Id such that their union equals {1,2,..., m} (or equivalently, such 
that l\ + I2 + ■ ■ ■ + Id = tti). The support supp(6) of a real vector b G W 71 is the subset of indices 
corresponding to the nonzero entries of b. The following condition requires that a certain linear 
subspace has a basis with disjoint supports. 

Condition 3.1. For a chemical reaction system given by a network G with m complexes and 
reaction rate constants k*-, let £ denote its complex-to- species rate matrix (|2.2p . and set d := 
dim(ker(£)). We say that the chemical reaction system satisfies Condition I3.il if there exists a 
partition I\, I2, ■ ■ ■ , Id of {1, 2, . . . , m} and a basis b 1 , b 2 , . . . , b d G W 71 of ker(S) with supp(6') = I{. 

Remark 3.2. Conditions 13.11 1 3.4| and 13.61 in this article are essentially linear algebra conditions. 
When we consider a specific choice of rate constants k*-, checking these conditions involves com- 
putations over M. However, the objects of interest (such as the subspace in Condition 13. If) are 
parametrized by the unknown rate constants Kij, so verifying the conditions can become quite 
complicated for large networks. In this case, we need to do linear computations over the field 
Q)(kij) of rational functions on these parameters and check semialgebraic conditions on the rate 
constants (cf. Remark |3.T() . 

Condition 13.11 implies that the steady state ideal Js* is binomial: 

Theorem 3.3. Consider a chemical reaction system with m complexes, and let d denote the di- 
mension o/ker(S). Assume that Condition \3. 1\ holds (i.e., there exists a partition ■ ■ ■ ,Id °f 
{1, 2, . . . , m} and a basis b 1 , b 2 , . . . , b d € M m of ker(S) with supp(6 4 ) = Ii). Then the steady state 
ideal Js* is generated by the binomials 

b^ x Vj2 - Jy> j2 x v n } for all ji, j 2 G Ij, and for all 1 < j < d. (3.1) 

Proof. Consider the vectors & - a = b J j i ej 2 — bj 2 eji G ^ m for all j%, J2 G Ij, for all 1 < j < d. It is 
straightforward to check that these vectors span the orthogonal complement ker(S)- 1 - of the kernel 
of S. But by definition, this complement is spanned by the rows of the matrix S. Therefore, the 
binomials 6^^2(2;) — b 3 ^j x (x) are M-linear combinations of the polynomials fi(x), f2(x), . . . , f s (x), 
and vice- versa. And so the binomials in (13. ip give another system of generators of ■ □ 

Note that Theorem 13.31 does not provide any information about the existence of (toric) steady 
states (i.e. real solutions to the binomials (|3.ip . cf. Definition 12. 2p . let alone positive steady states. 
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In general, this is a question of whether a parametrized family of polynomial systems has real 
solutions. For this purpose two further conditions are needed: 

Condition 3.4. Consider a chemical reaction system given by a network G with m complexes and 
reaction rate constants n*j that satisfies Condition 1 3. l\ for the partition 1%, I2, ■ ■ ■ , Id of {1, 2, . . . , m} 
and a basis b 1 ,b 2 , . . . ,b d G W 71 of ker(E) (with supp(6 l ) = U). We say that this chemical reaction 
system additionally satisfies Condition \3.4\ if for all j G {1,2, ... ,d}, the nonzero entries of V 
have the same sign, that is, if 

sign r&jj = sign f&Q , for all j 1} j 2 G Ij, for all 1 < j < d. (3.2) 
The next result can be used to check the validity of Condition [37 



Lemma 3.5. Consider a chemical reaction system with m complexes that satisfies Condition [TO] 
for the partition J 1; I2, ■ ■ ■ ,Id of {1,2, ... ,m} and the basis b l ,b 2 ,...,b d G W 71 of ker(E). Let 
j G {1,2,... ,d}, There exists an (Ij — 1) x Ij submatrix Ej of E with columns indexed by the 
elements of Ij and linearly independent rows (that is, rank(T,j) = Ij — 1). Let Ej be any such 
matrix. For i G {1, . . . ,lj}, call Ej(i) the submatrix of Ej obtained by deleting its i-th column. 
Then the system satisfies Condition \3.4\ (that is, equations (|3.2|) are satisfied) if and only if, for all 
j G {1, 2, . . . , d}, the sign o/det(Ej(i)) is different from the sign o/det(Ej(i + 1)) for 1 < i < Ij — 1. 

Proof. First, note that the kernel of the submatrix of E formed by the columns indexed by Ij has 
dimension one and is spanned by the vector b'j which consists of the Ij entries of V that are indexed 
by Ij. So there exist Ij — 1 rows that give a matrix Ej as in the statement. 

By a basic result from Linear Algebra, the kernel of Ej is spanned by the vector v' with i-th entry 
equal to (—1)* det(Ej(i)). As the vector b'j must be a multiple of v' , it is immediate that (|3 . 2[) holds 
if and only if the sign of det(Ej(i)) is different from the sign of det(E,-(i+ 1)) for 1 < i < L — 1. □ 

Condition 13.41 is necessary for the existence of positive real solutions to the system defined by 
setting the binomials (|3. If) to zero. In working towards sufficiency, observe that the system can be 
rewritten as 

x y n y h = -4i, for all j\, j 2 G Ij and for all 1 < j < d. 

bi 

n 

Note that Condition 13.41 implies that the right-hand side of the above equation is positive. In 
addition, we are interested in positive solutions x G M^q, so we now apply In (•) to both sides and 
examine the solvability of the resulting linear system: 

In x (yj 1 — yj 2 ) = In , for all j\ , j'2 G Ij and for all 1 < j < d, 
b n 

where Inx = (ln(xi), ln(x2), • • • ,ln(x s )). Now collect the differences (yj x — yj 2 ) 1 a s columns of a 
matrix 

A : = [(Vh ~ y^% 1)js el^i<j<d ' ( 3 - 3 ) 

and define the (row) vector 

e K :=(ln|) . (3.4) 

Observe that the basis vectors V and hence the vector Q K depend on the rate constants. The 
binomials (|3.ip admit a real positive solution (in the presence of Condition I3.4[) . if and only if the 
linear system 

(lnx)A = 6 K (3.5) 
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has a real solution (In a;) G M s . This is the motivation for our final condition and Theorem 13.81 
below: 

Condition 3.6. Consider a chemical reaction system given by a network G with m complexes and 
reaction rate constants k|- that satisfies both Condition \3.1\ (i.e. there exists a partition I±, I2, ■ ■ ■ , Id 
of {1,2,..., m} and a basis b 1 , b 2 , . . . , b d G W m of ker(S) with supp(6 J ) = U) and Condition \3.4\ 
(i.e., the coefficients of each binomial in equation h3. 1\) are of the same sign). Recall the matrix A 
and the vector K (defined in equations \3. 3\) and |5^P, respectively). Let U be a matrix with 
integer entries whose columns form a basis of the kernel of A, that is, U is an integer matrix of 
maximum column rank such that the following matrix product is a zero matrix with s rows: 

AU = . 

We say that this chemical reaction system additionally satisfies Condition \3.6\ if the linear sys- 
tem (|3.5p has a real solution (lnx) G W. Equivalently, the Fundamental Theorem of Linear Algebra 
[31J implies that equation (|3.5p has a solution, if and only if 

e K U = 0. (3.6) 

Remark 3.7. Conditions 13.41 and 13.61 impose semialgebraic constraints on the rate constants: 

• If the matrix A defined in (I3.3P has full column rank (i.e. the right kernel is trivial), then U 
is the zero vector. It follows that equation (]3,6p holds, and hence, Condition 13.61 is trivially 
satisfied for any positive vector of rate constants. We will see that this is the case for 
multisite phosphorylation networks. 

• If the matrix A does not have full column rank (i.e. there exists a nontrivial right kernel), 
then equation (j3.6j) can be translated to a system of polynomial equations in the rate 
constants. 

Now we can state sufficient conditions for a chemical reaction system to admit positive toric 
steady states: 

Theorem 3.8 (Existence of positive toric steady states). Consider a chemical reaction system 
with m complexes which satisfies Condition \3.1\ and hence has a binomial steady state ideal Jj^, . 
Then this chemical reaction system admits a positive toric steady state if and only if Conditions \3.4\ 
and \3.6\ hold. 

Proof. Assume that Conditions 13 . 1 (. [3.41 and 13.61 hold. Lemma 13.51 implies that the coefficients of 
the binomial system are of the same sign, hence A and Q K given in equations (]3.3p and (I3.4p and 
the linear system (|3.5p are well-defined. Then Condition 13.61 gives a solution (In x) to the system 
(|3.5p . which immediately gives a positive steady state x G K> of the chemical reaction system. 

On the other hand, assume that Condition 13 . 1 1 holds and that the system admits a positive steady 
state, that is, the binomial system (|3.ip has a positive real solution. In this case the coefficients of 
the binomials must be of the same sign, which implies that Condition 13.41 holds additionally. Again, 
positive real solutions of the binomial system imply solvability of the linear system (|3.5|) and thus, 
Condition 13.61 is satisfied as well. □ 

Remark 3.9 (Existence of steady states using fixed point arguments). In some cases, one can 
establish the existence of positive steady states by using fixed-point arguments. If the stoichiometric 
compatibility classes of a network are bounded, a version of the Brouwer fixed point theorem 
guarantees that a non-negative steady state exists in each compatibility class. If moreover the 
chemical reaction system has no boundary steady states, we deduce the existence of a positive 
steady state in each compatibility class. For example, the multisite phosphorylation networks 
that are studied in this article have this property. The positive conservation laws in (|4.ip ensure 
boundedness and Lemma [431 shows that no boundary steady states can occur. 
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The focus of our results, however, is slightly different. We are more interested in parametrizing 
the steady state locus (and hence all positive steady states) and less with the actual number of 
steady states within a given stoichiometric compatibility class (apart from Section [5j where we are 
concerned with compatibility classes having at least two distinct positive steady states). Moreover, 
using fixed point arguments, the existence of positive steady states may only be deduced if the 
chemical reaction system has no boundary steady states, which is somewhat rare in examples from 
Computational Biology. Our results do not require any information about boundary steady states. 

Example 3.10 (Triangle network, continued). We return to Example \2.3\ to illustrate the three 
conditions. First, ker(S) is the plane in M 3 orthogonal to the vector (— 2k\2 — ^13, 2^21 + K 23, K 3i — 
^32)- It follows that the partition {1,2}, {3} works to satisfy Condition \3.1\ if and only if K31 = K32. 
Therefore, for a chemical reaction system arising from the Triangle network, Condition \3.1\ holds 
(with partition {1, 2}, {3} ) if and only if the system has toric steady states. The forward direction 
is an application of Theorem \3.3l while for general networks the reverse implication is false: we will 
see in Example \3.15\ that there are networks with toric steady states that do not satisfy Condition \3.1\ 
for any partition. 

Next, for those systems for which K31 = ^32, Condition \ 3.4\ comes down to verifying that the 
entries of the vector (—2k\2 — K13, 2«2i + K23) have opposite signs, which is clearly true for positive 
rate constants. Finally, Condition \3.6\ asks (again, in the K31 = K32 setting) whether the following 
linear system has a real solution (lnxi, In £2) G M 2 : 

(lux,, lnx 2 )^j2 = ! n (StS) ' 

=A =0 K 
which is clearly true. This linear equation arises from the binomial equation 

(2«12 + K13) x\ - (2k 2 i + K23) x\ = 0. 
As Condition \3.6\ holds. Theorem 1 3. 8\ implies that these systems admit positive steady states. 

Under the hypothesis of Theorem 13.8} the following result shows how to parametrize the steady 
state locus. 

Theorem 3.11. Consider a chemical reaction system that satisfies Conditions \3.1\ \3.4\ and \3.b\ 

Let A G Z wxs be a matrix of maximal rank w such that ker(^4) equals the span of all the differences 
Vh ~ Vji f or ii> J2 G Ij> where 1 < j < d. For 1 < i < s, we let A4 denote the i-th column of A. 
Let x G M^q be a positive steady state of the chemical reaction system. Then all positive solutions 
x G M> t° th e binomial system )13. 1\) can be written as 

x = ( Xl t A \ x 2 t A \ x s t A °) , (3.7) 

for some t G M> (where we are using the standard notation for multinomial exponents). In 
particular, the positive steady state locus has dimension w and can be parametrized by monomials 
in the concentrations. Any two distinct positive steady states x 1 and x 2 satisfy 

lux 2 - lux 1 G im [A 1 ) = span {y j2 - y jl \ ji,j 2 G Ij, 1 < j < d}^ . (3.8) 

Proof. By definition, the rows of A span the orthogonal complement of the linear subspace spanned 
by the differences yj 2 — yj 1 for ji,j 2 G Ij, 1 < j < d. Let x G M> be a positive steady state of 
the chemical reaction system; in other words, it is a particular positive solution for the following 
system of equations: 

tr> x yj2 - br>.x Vj i = for all 71,72 G L, and for all 1 < j < d . 

(Here the b 3 are the basis vectors of ker(S) with disjoint support.) Then it follows from basic 
results on binomial equations that all positive solutions x G Rt, to the above system of binomial 
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equations can be written as 



x = {xxt Al , x 2 t M , x s t As ) , 



for some t € K>o- In particular, the positive steady state locus has w degrees of freedom. 
For the convenience of the reader, we expand now the previous argument. In fact, it is easy to 
check that any vector of this shape is a positive solution. We first let x* be a particular positive 
solution of the above binomials. Then 3=- := ( . . . ^f- ) is a positive solution of the system 

of equations: 

x yj 2 - x Vji = o for all jx,j 2 G Ij, for all 1 < j < d . 

Therefore, (^) v ^~ y n = i. Or, equivalently, m(C) . (y h - y h ) = 0. This implies that In [€) 
belongs to the rowspan of A, and this means there exist Ai, A2, • • • , A w such that, if Ai, A 2 , ■ ■ ■ , A 
represent the rows of A, then we can write 

ln(^-j) =Ai(^4i)i + A2(-42)i + --- + A 1JJ (A )i, for all 1 < i < s . 



w 



X 

If we call ti := exp(A^) for 1 < £ < w, then x* = X{t Ai for all 1 < i < s, which is what we wanted 
to prove. □ 

We now turn to the case of a network for which Condition 13.11 holds with the same partition for 
all choices of rate constants. The following result, which follows immediately from Theorem 13.81 
states that for such a network, the semialgebraic set of rate constants that give rise to systems 
admitting positive steady states is defined by Conditions 13.41 and 13.61 

Corollary 3.12. Let G be a chemical reaction network with m complexes and r reactions, and 
assume that there exists a partition I\ , I 2 ■ ■ ■ , Id of the m complexes such that for any choice of 
reaction rate constants, the resulting chemical reaction system satisfies Condition \3.1\ with this 
partition. Then a vector of reaction rate constants k*- £ M'^q gives rise to a system that admits a 
positive steady state if and only if K*j satisfies Conditions \3.4\ and \3.6[ 

In the following example, we see that the 2-site phosphorylation network satisfies the hypothesis 
of Corollary 13.121 The 2-site system generalizes the 1-site system in Example 12. 1\ and we will 
consider general re-site systems in Section [H 

Example 3.13 (2-site phosphorylation system). The dual phosphorylation network arises from the 
1-site network (|2.4p by allowing a total of two phosphate groups to be added to the substrate of So 
rather than only one. Again there are two enzymes (E and F), but now there are 3 substrates (So, 
Si, and S 2 ). The substrate Si is the substrate obtained from So by attaching i phosphate groups 
to it. Each substrate can accept (via an enzymatic reaction involving E) or lose (via a reaction 
involving F) at most one phosphate; this means that the mechanism is "distributive" . In addition, 
we say that the phosphorylation is "sequential" because multiple phosphate groups must be added in 
a specific order, and removed in a specific order as well. 



S + E^ ES c 4° Si+E^ ES X °4 l S 2 + E 

^offo ^offi 

(3.9) 

^onj 1 ^onn , 

S 2 + F ^ FS 2 C 4 X Si + F i=* FS 1 c 4° So + F 

loff 1 ^off 

We order the 9 species as (So, Si, S 2 , ESq, ESi, FS±,FS 2 , E, F), and we order the 10 complexes as 
(So + E,S! + E, S 2 + E, ESq, ESi, S + F, Si + F, S 2 + F, FS U FS 2 ). The 9 x W-matrix Y* and 
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the 10 x 10-matrix A l R for this system are the following: 
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We will analyze the steady state locus of the resulting chemical reaction system by focusing on the 
structure of the kernel of the matrix £ = Y t A t K of the system. Note that the network 113. 9\) has only 
two terminal strong linkage classes, {S2 + E} and {So + F}. Also, span{e3,eg} C ker(E), where ej 
denotes the i-th canonical vector 0/R 10 . A partition of the 10 complexes that satisfies Condition \3.1\ 
is given by I\ = {1,4,7,9}, I2 = {2,5,8,10}, I3 = {3}, and I4 = {6}. A corresponding basis of 
ker(E), that is, one in which the i-th basis vector has support Ii, is: 



b 1 



off 



^cato )^oni ^cati ^oni ^ono ^cato 




^oriQ ^oni k ca ± 1 ^oni ^ono ^cato 






^ono ^catg ^0111 ^cati ^oni (^cato 



kff ) 



^-onQ ^cato ^ono ^oni ^cati ^oni 





^ono ^cato ^ono (^off 1 H~ ^cati )^oni ^cati 




^oiio ^catg ^ono ^oni ^oni fcati 




^ono ^cato ^ono ^oni ^cati (^cati H~ ^offi) 



^ono ^cato ^ong ^oni &cati ^oni 



e3 



e6 



10 ; v 6 ker(E) if and only if v satisfies the 

(3.10) 



The structure of this basis {b 1 } implies that for v 6 
following binomial equations: 

b\v4 - b\v\ = , hjv 5 - b\v 2 = , 
b\v 7 - b\v x = , b%v 8 - bjv 2 = , 
b\v 9 - b\v\ = , bjvio ~ bl v 2 = , 

Hence, any steady state of the 2- site phosphorylation system must satisfy the following equations 
in the species concentrations x = (xs , a^si, . . . , xe, xf): 

b\x4 — b\xsXi = , &2 X 5 ~~ t>l X 8 X 2 = ® ' 

b\xQX 2 — bjx$x± = , b\xo,x% — b\x)i,X2 = , (3-11) 



b\xQ — bgXgXi = , &2 X 7 — b1 Q X8X2 = . 

To check Condition \3.(j\ we consider the matrix A and the vector Q K : 
A = [e 4 - e 8 - ei | e 9 + e 2 - e 8 - e\ | e 6 - e 8 - ei | e 5 - e 8 - e 2 | e 9 + e 3 



e 2 e 7 



e-2\ 
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/, b\ , bl , bl , bl , 6g , o? n 
B = In — In — In — In — In — In — 

K \ b\' b\ : b\' &1' bl 

It is straightforward to check that A has rank 6 and hence full rank. Thus Condition \3.6\ is trivially 
satisfied and does not pose any constraints on the rate constants. 

Following the proof of Theorem \3.11l we first will parametrize the solution set of the following 
reduced system: 

x 4 - x 8 x 1 = , x 5 - x 8 x 2 = , 

X9X2 — x$x\ = , X9X3 — x 8 X2 = , (3-12) 

Xq — X 8 X\ = , £7 — X 8 X2 = . 

We are interested in an integer matrix A such that kei(A) = im (A). One such matrix is 

12 12 12 10 

.1 1000111111 

11111110 

This provides the following 3- dimensional parametrization of the reduced system: 

(tl,t2,ts) !->■ (*3, tits, tft3, M2^3) £1*2*3, *1*2*3? *1*2*3) *1*2; h) , 

where t2 is the concentration of the enzyme F, ti is the quotient of the concentration of the enzyme 
E divided by the concentration of the enzyme F, and t% is the concentration of the substrate Sq. 
Returning to the original binomials (|3.1ip . we have the following particular solution: 

x 1 -x 8 -x 9 -i, x 2 -^, s 3 - 6 i 6 2, a;4- 6 i» ^5- H - ^ , x 7 - ^ . 

Therefore we obtain the following 3- dimensional parametrization of the positive steady state locus 
of (|3.1ip . as predicted in Theorem \3.11\ 

M 3 >0 ^M 9 >0 (3.13) 
, ( b\ b 2 8 b\ 2 b\ b\b\ 2 b\ b\j>\ 2 \ 

(il,*2)*3j ^ I *3) ^T*l*3) ^2*1*3) J^hhhi ^2*1*2*3, p-*l*2*3, ^1^2 n*2*3; *1*2, *2 I • 
\ 1 12 1 12 1 12 ^ 

Recall that the values 6*- are polynomials in the rate constants shown in the display of the vectors 
b 1 and b 2 . Finally, note that none of the calculations in this example depends on the specific 
values of the rate constants; in particular, one partition works for all systems, so the hypothesis of 
Corollary \3.12\ holds. 



3.1. More general sufficient conditions. We show in Example 13. 151 below, extracted from [30J , 
that the conditions in Theorem 13.81 are not necessary for a chemical reaction system to have toric 
steady states; in other words, the converse of Theorem 13.81 does not hold. However, the condition 
for the steady state ideal to be generated by binomials always can be checked algorithmically via 
a Grobner basis computation, as stated in the following lemma. 

Lemma 3.14 (Proposition 1.1. (a) of [ID]). Let I be a binomial ideal, let -< be a monomial order, 
and let G be the reduced Grobner basis of I for that ordering. Then G consists of binomials. 

Lemma 13.141 is a basic result about binomial ideals which is due to Eisenbud and Sturmfels |10j ; 
it is a result concerning polynomial linear combinations. Note however that Theorem 13.81 requires 
only linear algebra computations over K. We make use of Lemma 13.141 in the following example. 
We will return to it later to show that Theorem 13.191 below can be used to prove that this system 
has toric steady states, without needing to compute a Grobner basis. 
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Example 3.15 (Shinar and Feinberg network). This example demonstrates that Condition \3.1\ 
is not necessary for a chemical reaction system to have toric steady states. The network in Ex- 
ample (S60) of the Supporting Online Material of the recent article of Shinar and Feinberg is the 
following [30j: 

XD +± X +± XT X p 
X p + Y <=> X P Y K -% X + Yp 
XT + Y p £± XTYp K ^>° XT + Y ( 3 - 14 ) 

K11 ' 12 «,,„ 

XD + Y p +± XDYp 4 XD + Y 

«12,11 

We denote by x\, x 2 , ■ ■ ■ ,Xg the concentrations of the species as follows: 

XXD = Xl, X X = X 2 , X X T = X 3 , X Xp = X 4 , 
XY = X5, XX P Y = XQ, X Yp = X 7 , XXTY P = Xg, X X DY P = Xg . 

Note that the numbering of the 13 complexes in the network is reflected in the names of the rate 
constants Kij . The chemical reaction system is the following: 

= -K12X! + K21X2 ~ K\\,\2X\X 7 + («12,11 + «12,13)a:9 

-p. = Kl2 x x + (-K21 - K 2 ?,)X2 + ^32^3 + ^67^6 

~W = K 23 X 2 + (-K32 - K U )x 3 - KggX 3 X 7 + (K9S + Kg iW )x 8 

= Ku x 3 - K 56 X 4 X 5 + K 65 X 6 

% = -^56^4^5 + ^65X6 + ^9,10^8 + ^12,13^9 (3-15) 

■4^ = K 56 X 4 X 5 + (-«65 - ^67)^6 

% = K 6 7^6 - ^9X3X7 + K 98 X 8 - Kn tl 2XiX 7 + K 12yU Xg 

= KggX 3 X 7 + (-«9 8 - K 9i i )x 8 

-jf = Kn,12XlX 7 + (-Kl2,ll - «12,13)iC9 

The reduced Grobner basis with respect to the lexicographical order x\ > x 2 > X4 > X5 > X6 > x 8 > 
xg > X3 > X7 consists of the following binomials: 



91 = 


[k89K12K23K 9j io(k12,11 + «12,13) + Kll,12K2lKl2, 13(^98 + K 9i i )(k32 + K34)]x 3 0;7 + 
+ [-«23«34«12(«12,ll + «12,13)(«98 + Kg >W )]x 3 


92 = 


[-K 11 .i 2 K 2 iK 3 4{k 9s + K 9j i )(k32 + K34)]x 3 + 

+ [«11,12«21«12,13(«98 + K9,io)(«32 + K34) + Ki 2 K 23 K 89 Kg, io(«12, 11 + Kl2,13)]^9 


93 = 


[-K23K34K 89 K12(ki2,11 + K 12 ,13)]x 3 + 

+ [K23K 9 ,10K89Kl2(Kl2,ll + «12,13) + Kll,12K21«12, 13(«98 + K9,w)(k 32 + K 3 4)]x S 


94 = 


K67X6 - H 3i X 3 


95 = 


K56H67XiX 5 + K 3i (-KG5 ~ Kq 7 )x 3 


96 = 


K 23 X 2 + (-K32 - K34)x 3 


97 = 


-k 2 i(k 32 + k 34: )x 3 + H 12 K 23 Xx 



(3.16) 



Therefore, the network has toric steady states (for any choice of positive reaction rate constants ) 
because the steady state ideal can be generated by g\, g 2 , g 7 . However, we claim that this 

chemical reaction system does not satisfy Condition \3.1[ In fact, for any rate constants, it is not 
possible to find a partition I\, I 2 , ■ ■ ■ , Ie Q {1, 2, . . . , 13} such that ker(E) has a basis {b l ,b 2 , . . . , ft 6 } 
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with supp(6 J ) = 1{. This can be seen by noting that the kernel o/S can be generated as follows: 

, /v,N / ( «2lKl2. 13(^32 + K U )\ ( K12 ,13(^32 + «34)\ . ( Kl2,is\ . , Q i-s 

ker(S) = ( e 4 , e 7 , ei , ei 3 , — e i + e 2 + — e 3 + (3.17) 



K23K34K12 / \ K23K34 / \ K34 

(k 6 5 + K67)K12.13\ , />«;L2,13\ , f (Kl2,ll + «12,13) A 

e 5 + e 6 + en + ei2, 



^67^56 / \ K 67 / \ ^11,12 

K2lKg,io(K32 + K 34 )\ . / Kg i (k 32 + K34) \ , /K9,lo\ 

ei + e 2 + e 3 - 



K 23 K 34 K 12 / V K 23 K 34 / \ K 34 

/ («65 + K67)K9.1o\ . /«9,lo\ . /«98 + «9,lo\ 

+ e 5 + — e 6 + e 8 + e 9 

V K67K56 / V «67 / V K 89 / 

Our next result, Theorem 13.191 will generalize Theorem 13.81 by giving a stronger condition that 
guarantees that the steady state locus is generated by binomials. We first need to generalize 
Conditions 13.11 13.41 and 13.61 to any (finite) polynomial system. 

First we must introduce some notation. For polynomials F±, F2, . . . , F s i G M.[xi,X2, ■ ■ ■ ,x s ], we 
denote by x Vl , x V2 , . . . , x Vm ' the monomials that occur in these polynomials; that is, there exist 
F{j G R such that F{(x) = Y^j=x ^%j xVl for i = 1,2, s'. We can write the polynomial system 
F x (x) = F 2 (x) = ■■■ = F s ,{x) = as 

T! ■ V'(x) = , (3.18) 
where S' = (Fy) G R s ' xm ' i s the coefficient matrix and W(x) = {x Vl ,x V2 ,. . . ,x v ™') 1 . We will let d! 
denote the dimension of ker(E'). 

Condition 3.16. We say that the polynomial system (|3.18|) satisfies Condition \3.16\ if there exists 
a partition I\, I2, ■ ■ ■ , Id' of {1, 2, . . . , m'} and a basis b 1 , b 2 , . . . , b d G R m of ker(S') such that 
supp(6 l ) = Ii. 

Condition 3.17. Consider a polynomial system (|3.18p that satisfies Condition \3.16\ for the parti- 
tion h, h, Id' of {1,2, ... ,m'} and a basis b 1 , b 2 , . . . , b d ' G W m ' o/ker(S') (with supp(6 l ) = I{). 
We say that the system satisfies additionally Condition \3.ll\ if for all j G {1,2, ... ,d'}, the 
nonzero entries of b 3 have the same sign. 

As before, we collect the differences of exponent vectors as columns of a matrix 

A' := [(y h - yn)% n , j2 ei 3 ,yi<j<d' ( 3 - 19 ) 

and define the (row) vector 

8':=(ln^) . (3.20) 

Condition 3.18. Consider a polynomial system (|3.18p which satisfies Conditions \3.16\ and 3.11 . 
Let V be a matrix with integer entries whose columns form a basis of the kernel of A'. We say 
that this system satisfies additionally Condition \3.18l if the following holds: 

@'U' = . 

We then have the following sufficient conditions: 

Theorem 3.19. Consider a chemical reaction system with m complexes and assume that there 
exist monomials x ai ,x a2 , . . . , x ae and indices i\, ii, ■ ■ ■ , it, with {i\, ii, ■ ■ ■ , ie} Q {1,2, ... ,s}, such 
that Condition \3.16\ holds for the enlarged polynomial system 

f 1 = ... = f s = x ^f il = ... = x ^f k =0. 

Then the steady state ideal is binomial. 

Moreover, the system has positive (toric) steady states if and only if Conditions 3.11 and \3.l8\ 
hold additionally for the enlarged system. 
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This theorem can be proved following the lines of the proof of Theorem 13.81 for the enlarged 
system defined in the statement. It is important to note that the ideal (fi, f 2 , . . . , / s ) equals the 
ideal </!,.. .,f s ,^f h ,...,^f k ). 

With similar proof as in Theorem I3.11| we moreover have: 



Theorem 3.20. Under the hypotheses of Theorem \3.19l the steady state locus can be parametrized 
by monomials in the concentrations. 

Remark 3.21. As with Conditions 13.11 [ 3.4| and 13. 6\ checking the Conditions in the statement of 
Theorem 13.191 involves linear algebra computations over R for fixed rate constants or over Q(kij) 
for generic parameters, once the monomials x ai are given. In small cases, such monomials can be 
guessed. In the following example, they were traced in the standard algorithm for the computation 
of a Grobner basis of the ideal (f%, . . . , f s ). 



We end this section by returning to Example 13.151 
Example 3.22 (Shinar and Feinberg network, continued). Consider the system of equations: 



( h 

h 








/ 9 = o 

xjfi = 
x 7 h = 

X7fs = 

I x 7 f 9 = 



(3.21) 



This enlarged system satisfies Conditions 1 3. 16\ and 3.17\ for the following partition: 

h = {4}, h = {10}, h = {13}, h = {14, 15}, I 5 = {16, 17}, h = {1, 2, 3, 5, 6, 7, 8, 9, 11, 12} 
and the following basis b , 6 2 , . . . , 6 6 of its kernel verifying supp(6^ 



If 



b 1 =e 4 , b 2 = do , b 
b e 



ei3 , b 4 = (fci2,n + fci2,i3)ei4 + fcn,i 2 ei5, b 5 = (fc 98 + fcgio)ei 6 + fc 89 ei 7 , 

(^12^23^89^9, lo(&12,ll + fel2,13) + ^21^11,12^12,13(^32 + k 34 )(k 98 + k 9 . w ))k 2 l (fc 32 + k 3i )k 56 k 67 ei + 
(^12^23^89^9, 10(^12, 11 + ^12,13) + fel^ll, 12^12, 13(^32 + ^34X^98 + &9,lo))&12 (k 3 2 + k 34 )k 56 k 67 e 2 + 
(^12^23^89^9, 10(^12, 11 + &12,13) + ^21^11, 12^12, 13(^32 + k 3i )(k 9 $ + k 9 ^ ))ki 2 k2 3 k 56 k 67 e 3 + 
(^12^23^89^9, 10(^12, 11 + fcl2,13) + k 2 ikn, 12^12, 13^32 + k 3i )(k 9s + k 9A0 ))ki 2 k 23 k 34 (k 65 + ^67)65 + 
(^12^23^89^9, 10(^12, 11 + &12,13) + fel^ll, 12^12, 13(^32 + fc34)(^98 + k 9 ^ ))ki 2 k 23 k 3 4k 5 6e 6 + 
k\ 2 k 23 k 3i (k 32 + k 34 )k 56 k 67 (k 9s + fcg,io)(&12,ll + ^12,13)67 + 

^12^23^34^56^67(^98 + &9,lo) (^12,11 + &12,13)es + ^12^23^34^56^67^89(^12, 11 + &12,13)e9 + 

k\ 2 k 2 ik 23 k 3i {k 32 + k 3i )k 56 k 67 (k 9 $ + fcg,io)(fci2,ii + fci2,i3)en + 

k\ 2 k 2 ik 23 k 3 i{k 32 + ^34) &56 &67 (&98 + fc9,lo)fcll,12ei2 • 



In addition to the monomials already occurring in f±, f 2 , . . . , /g, the following 4 monomials are 
also in the augmented system: x Vl4 = X\x 7 , x Vl5 = Xgx 7 , x Vl6 = X3X7, and x yi7 = x$x 7 . By 
Theorem \3.19l the system has toric steady states. Recall that the binomials gi, g2, ■ ■ ■ , g 7 in equa- 
tion (|3.16p generate the ideal (fx, f 2 , ...,/ 9 ) = f 2 , • 2:7/1, x 7 f 3 , x 7 f 8 , x 7 f 9 ). We can 
see immediately that there are positive steady states for any choice of positive rate constants, and 
so there is no need to check Condition \3.18[ 
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4. THE n-SITE PHOSPHORYLATION SYSTEM HAS TORIC STEADY STATES 

In this section we introduce the n-site phosphorylation system (under the assumption of a dis- 
tributive and sequential mechanism). To show that these systems have toric steady states, we 
apply Theorem 13.8) this generalizes Example 13.131 (the n = 2 case). Further, we note that the 
parametrization of the steady state locus given by Theorem 13.111 is implicit in work of Wang and 
Sontag [36 . 



4.1. The n-site phosphorylation system. We now define the n-site phosphorylation system 
(also called a "multiple futile cycle") S n (K,C), which depends on a choice of rate constants k £ R>*q 
and values of the conservation relations C = (-Etot> ^tot> ^tot) e ^>o- As i n * ne earlier example of 
the 1-site network (12. 4h and the 2-site network (13. 9p , we will make the assumption of a "distributive" 
and "sequential" mechanism (see, for example, [6]). As discussed in the Introduction, this n-site 
phosphorylation system is of great biochemical importance: it is a recurring network motif in many 
networks describing processes as diverse as intracellular signaling (e.g. MAPK signaling with n = 2 
and n = 3), cell cycle control (e.g. Sicl with n = 9), and cellular differentiation (e.g. NFAT with 
n = 13). 

Following notation of Wang and Sontag [36], the n-site phosphorylation system arises from the 
following reaction network: 



So + Ep ES Sl + E S\ + F ™ F Sl H° S + F 

fc cat n _! fonn-i W 



S n + F ^ FS n -T 1 S n _i + F 



We see that the n-site network has 3n+3 chemical species So, . . . , S n , ESo, . . . , ES n -\, FSi, . . . , F5 n , 
F, and F, so we write a concentration vector as x = (sq, . . . , s n , Co, . . . , c n _i, di, . . . , d n , e, /), which 
is a positive vector of length 3n + 3. These species comprise 4n + 2 complexes, and there are 6n 
reactions. Each reaction has a reaction rate, and we collect these in the vector of rate constants 
K = (k ono , . . . , Zcat n _i) € K^Jq. 

For our purposes, we will introduce the following numbering for the complexes (which is com- 
patible with the numbering in Examples 12.11 and 13. 13[) : 



1 ^ n+ 2 -> 2 

2 n + 3 -> 3 



n<^2n + l— >-n+l 



2n + 3 ^ 3n + 3 2n + 2 
2n + 4 +± 3n + 4 -> 2n + 3 



3n + 2 ^ 4n + 2 3n + 1 



The conservation relations here correspond to the fact that the total amounts of free and 
bound enzyme or substrate remain constant. That is, the following three conservation values 
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C = (-^tot > -^tot > ^tot ) G K^q remain unchanged as the dynamical system progresses: 

n-l 

^tot = e + Y Ci ' 

i=0 
71 

*tot = /+!>> f 4 - 1 ) 

i=i 

n n—1 n 

s tot = Y s < ■ + Y Ci + Y di ■ 

i=0 i=o i=l 

Any choice of these three values defines a bounded stoichiometric compatibility class of dimension 
3n: 

Vc = |x G M^, +3 | the conservation equations (|4.ip holdj . 

Note that the right hand side of each of the three conservation relations (|4.ip is defined by a vector 
Zi € , i = 1, 2, 3. These vectors play an important role in the following remark and in Lemma [4, 2 1 
below. 

Remark 4.1 (Positive steady states by fixed point arguments). As indicated in Remark l3.9l one may 
deduce the existence of at least one positive steady state in each stoichiometric compatibility class 
Vc (defined by positive C) by fixed-point arguments, provided (i) Vc is bounded and (ii) there are 
no boundary steady states in any stoichiometric compatibility class Vc (with C G R>o)- Point (i) 
follows from the definition of Vc G M> given above. With respect to (ii), we point to Lemma 14.21 
below (which can be established by a straightforward generalization of the analysis due to Angeli, 
De Leenheer, and Sontag in Examples 1 and 2 in [TJ § 6]). 

Lemma 4.2. Let x* G M> — M> ^ e a boundary steady state. Set A := {i G {1, . . . , s} : x* = 0}. 
Then, A contains the support of at least one of the vectors G defining the conservation 
relations (|4.ip . Therefore, there are no boundary steady states in any stoichiometric compatibility 
class Vc with C G B^ - 

We will see in Theorem 14.31 that the steady state locus in this system is 3-dimensional. A 
forthcoming work will concern the question of how many times the steady state locus intersects the 
relative interior of a compatibility class Vc for multisite phosphorylation systems [3]. 



4.2. Results. For the n-site phosphorylation system, we will call its complex-to-species rate ma- 
trix E n , and we will let G n denote the underlying digraph of the chemical reaction network. In 
order to apply the results of Section [3] to this system, we now aim to exhibit a specific basis of 
the kernel of S n that satisfies Condition 13.11 We begin by describing the rows of S n := Y t ■ A* K 
as linear combinations of the rows of A K . Recall that A K is minus the Laplacian matrix of the 
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associated digraph. Letting R(i) represent the i-th row of A* K , we have: 

R(l) + R(2n + 2) 
R{2) + R(2n + 3) 

R(n + 1) + R(3n + 2) 



Y ■ A: 



R(ni 


-2) 


R(2n- 


M) 


R{3n- 


r3) 


i?(4n- 


r2) 



3n+3)x(4n+2) 



(4.2) 



+ #(2) + • • • + R(n + 1) 
R(2n + 2) + ■ ■ ■ + R(3n + 2) 

Our next aim is to exhibit a submatrix of S n that shares the same kernel as S n . The only relations 
that exist among the rows of A l K arise from the fact that the sum of the rows in each of the four 
blocks equals zero. Consequently, it is straightforward to check that 

rank(S ra ) = 3n . 

Moreover, if we delete any of the first 3n + 1 rows and the last two rows of S n , we obtain a new 
matrix that has maximal rank. As we are interested in describing the kernel of T, n , we will discard 
the first and the last two rows, and we will focus on the resulting submatrix. Furthermore, as the 
(n + l)-st and (2n + 2)-nd columns on S n are equal to zero, we already know that e n+ \ and e2n+2, 
the (n + l)-st and (2n + 2)-nd canonical basis vectors of M 4n+2 , belong to ker(S n ). Hence we can 
now focus on an even smaller submatrix of S n obtained by deleting the first and the last two rows, 
and the (n + l)-st and (2n + 2)-nd columns. We will call this submatrix T,' n , and we will denote by 
C(j) the column of T,' n which corresponds to the j-th column of S n after deleting the first row and 
the last two (for example, C(n + 2) will represent the (n + l)-st column of T,' n ). Then, if we call 
the submatrix of Y! n formed by its first 3n columns, the system T,' n v = is equivalent to the 
following one: 



t'3n 



C(3n + 3) 



C(4n + 2) 



V3n+1 



(4.3) 



Let us call 

D := det(S^) . (4.4) 
If D ^ 0, then we can use Cramer's rule to solve system (14.3p , In fact, we will show in Proposition l4.4l 
that this is the case and that we can find solutions to the system X n u> = such that all the nonzero 
entries have the same sign. 

Next we introduce a partition and a set of basis vectors b l that will be used to show that the n-site 
system satisfies Condition 13.11 The partition Ji, I2, . . . , I n +2 of {1, 2, . . . , 4n + 2} is the following: 



I j = {j, n + j + 1, 2n + j + 2, 3n + j + 2} (fori <j <n), I n+1 = {n + l}, I n+2 = {2n + 2}. (4.5) 



The entries in our vectors b l will be certain determinants. More precisely, let D^j) be minus 
the determinant of the matrix obtained by replacing C(£(j)) by C(3n + j + 2) in for £(j) = 
j, n + j + 1, 2n + j + 2, where 1 < j < n: 



D t (j) = - det 



C(l)\...\C(3n + j 



|C(3n + 2) 



(4.6) 
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Note that D, Dj, D n+ j + \, and D2 n +j+2, for 1 < j < n, define polynomial functions of k on ^ >0 - 
We will show in Proposition 14.41 that these functions D, Dj, D n+ j + i, and D 2n +j+2 are nonzero and 
have the same sign, for 1 < j < n. 

Now we may define the vectors b 1 ,b 2 , . . . ,b n of 1R^ +2 by: 



for 1 < i < 4n + 2, where 1 < j < n. 



' D 3 


if i = j 


D n +j+i 


ifi = n + j + l 


< D2n+j+2 


if i = 2n + j + 2 


D 


if i = 3n + j + 2 




V 


otherwise 


main result 


in this section. 



(4.7) 



Theorem 4.3. T/ie n-site phosphorylation system has toric steady states. The steady state locus 
has dimension 3 and can be parametrized by 



R 3 R 3 "+ 3 

/, , , \ . I . Ajn+3. , D 2n +3 D 3n+2 D n+2 D n+2 D 2n+ i 

(11,12,13) ^ i3, — p; — tits, — — — ... — — — t 1 t 3 , — — EiC2E3, — p; — ■•■ — p; — Cit2i3i 

D D JD2ri4-3 Dp. n + \ « 

— — tll2C3, •••! "p; p; ■ ■ ■ ~F> t \ t 2 l 3, tit 2 , 12 

U\ L) n L>i U n -x 

Moreover, the system satisfies Condition \3.1\ with the partition I\, I 2 , ■ ■ ■ , I n +2 described in (14, 5p 
and the basis {b 1 , . . . ,b n } U {e n +i> e 2n+2} where the vectors W are defined in (|4.7p and e n +\ and 
&2n+2 are the (n + l)-st and (2n + 2)-nd vectors of the canonical basis o/R 4n+2 . In addition, it 
satisfies Conditions \3.4\ and \3.6l 

In particular, 

-P2w+3 T>2n+3 -Pgw+j Dn+2 Dn+2 £ ) 2ra+l D D D2n+3 ^3k+1 , , 

' D t ' D 1 D„ ' "dT' £>!"'£)„' Di' ~D~ n D x " ' " D„_ x ' ' 

is an explicit positive steady state of the system. 

We remark that the parametrization given in the statement of this theorem, which is one of the 
possible parametrizations provided by Theorem I3.11|, gives systematically what Wang and Sontag 
obtained "by hand" in [36]. We note that the fact that this variety (the steady state locus) has a 
rational parametrization is a special case of a rational parametrization theorem for general multisite 
post-translational modification systems as analyzed by Thomson and Gunawardena 



4.3. Proof of Theorem 14.31 We start with the following proposition: 

Proposition 4.4. Let D be the determinant defined in (|4.4|) . and let Dj, D n +j+i, and L>2n+i+2 be 
as in (14.61) . for 1 < j < n. Then each polynomial function D, Dj, D n +j+\, D% n +j+2 '■ M^q ~~ ^ ^ f or 
1 < j < n, never vanishes, and these functions all have the same constant sign on R^j. 

Proof. For this proof, we will denote by R{i) the i-th row of the matrix obtained from A f K after 
deleting columns n + 1 and 2n + 2. (Note that this notation differs slightly from that introduced 
in equation (|4.2p .) The proof has two steps: first we demonstrate that D 7^ on the positive 
orthant, and then we show that the other functions Dj, D n +j+i, and Dm+j+2 are also nonzero on 
the positive orthant and that their signs coincide with that of D. 
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To prove that D ^ on R^g, we proceed by induction on n. First, if n — 1, we have: 



S'/ 



v on 





^cato 
"&off ~~ k c 

o 



<-on 




In this case, D = —k ono k ca ,t l ono ^ 0, as we wanted. 

For the n > 1 case, we suppose now that the D/0 result is valid for G n -i, the network of 
the (n — l)-site phosphorylation system. In order to visualize the calculations, we will reorder 
the rows and columns of E", placing C(l), C(n + 2), and C(2n + 3) as the leftmost columns, and 
R{2) + R(2n + 3), R(n + 2), and R(3n + 3) as the uppermost rows. We notice that this ordering 
does not alter the sign of the determinants, hence we can write 



D = det 



/ 





^cato 


^ono 




\ 




^ono 


— ^off ~~ ^cato 


















^ono 







V 











B 


) 



— /Cono^-cato^ono det (5) 



(4.8) 



where B is a (3n — 3) x (3n — 3)-submatrix of E^. This matrix B does not include either C(l),C(n + 
2), C(2n + 3), nor the first (n + l)-st or (2n + l)-st rows of E^. We next will see how the matrix B 
can be interpreted as the 3(n — 1) x 3(n — l)-matrix E^_ l5 the corresponding matrix of the smaller 
network G n -\. This interpretation will allow us to conclude by the inductive hypothesis that D / 
in the positive orthant. 

For the purpose of interpreting this submatrix of E" as the matrix of G n _i, it is important to 
note that the deletion of C(l),C(n + 2), and C(2n + 3) from is equivalent to calculating E'' 
after having deleted these columns from A f K before calculating E n . In turn, it is also equivalent 
to having deleted all the reactions that begin at the first, (n + 2)-nd and (2n + 3)-rd complexes 
of the network. Once we have additionally deleted the first, (n + l)-st, and (2n + l)-st rows (i.e. 
R(2) + R(2n + 3), R(n + 2), and R(3n + 3)), we obtain a new submatrix of E n whose entries we 
can rename as follows: 



fconj — : ^onj_!i ^offj — : ^oflF j _ 1 ) ^catj — : ^ C atj_ii 'onj — : ^on j _ 1 ; fcffj — : ^oflF j _ 1 ' ^catj — : ^catj_! ■ 

In fact, this new matrix is the corresponding complex-to-species rate matrix E^_ x for the network 
G n -±, with corresponding rate constants indicated by primes. We can also establish a correspon- 
dence between the nodes of the two networks: letting j' denote the j-th node of G n -±, then f 
corresponds to the following node of G n : 

j + 1 if 1 < f < n (complexes So + E, . . . , S n -\ + E in G n _i) 

j + 2 if n + 1 < j' < 2n (complexes ESq, . . . , ES n -2 in G n -\) 

j + 3 if 2n + 1 < / < 3n - 1 (complexes So + F, 5 n _i + F in G n -i) 

j + 4 if 3n < j' <4n-2 (complexes FS , ■ ■ ■ , FS n -i in G n _i) . 

From this correspondence, it follows that det(-B) equals det(E^_ 1 ), which is nonzero by inductive 
hypothesis, and therefore D ^ 0, which we wanted to prove. 

We now complete the proof by verifying the following claim: the polynomial functions Dj, 
D n+ j + i, Z?2n+j+2 never vanish, and they all have the same constant sign as that of D on R^, (for 
1 < j < n). 



j' corresponds to < 
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We first prove this claim for the case j = 1. We again reorder the entries of the matrices as 
described above, and as this ordering does not alter the sign of the determinants, we can write: 



D 



n+2 



D 



2n+3 



det 



det 



det 



/ 


lo&o 


^catg 


lon 




\ 







^offo — ^cat 












^cato ^offo 












V 











B 


) 



— — (fcoffo + fccat )ion icato det(-B) , 



/ 





l O S 






\ 






















^catp ^offp 









V 











B 


J 



/ 





^cato 


l o S 






h 


&off ^cat 
















^cato ^offo 





V 











B 



^■oiiq ^cato (^cato ^off'o 



det(B) , 



~S'n-i- As we already 



where B is the same matrix we described in equation (|4.8p . That is, B 
know that D ^ 0, we deduce that det(-B) ^ 0. By examining equation (|4.8p and the display above, 
we conclude that the claim is true for j = 1. 

For the j > 1 case, we will prove our claim by induction on n. The base case is n = 2 (as 
j > 1 is not possible when n = 1). In this case, the functions of interest are the following 
positive functions on Mi 2 : D 



^ong ^cato ^ono ^oni &cati ^oni > B2 



^ono^cato^ono (^offi ^cati )^oni ^cati j 



— ^ono ^cato^ono ^oni ^oni ^cati ; and D$ 



&onofccato^on fconi&cati(kati +^offi)- Hence our claim holds 



for n = 2. 

We now assume that the claim is true for G n -\. 



As we did above, we view G n -\ as a subgraph 



of G n , and if we call D' p( .,^ the corresponding determinant of the (n — l)-site system (for £(j') 



, (n - 1) + j' + 1, 2(n - 1) + j' + 2, for 1 < j' < n - l), then we have: 



for £(j')=i', 2(n- 

V)' 



^on ^cato^ono-D£(j') , (4-9) 

1) + j' + 2, where 1 < j' < n — 1. By the inductive hypothesis, the 
claim holds for the -D'^-n, so by equation (|4.9H . the claim holds for the -D^Q) as well. This completes 

□ 



the proof. 

We now take care of the zero entries of the vectors W defined in (|4.7p . We start by defining 
D u4rJfV as minus the determinant of the matrix obtained by replacing column C(u) by C{v) in 
for 1 < u < 3n + 2 such that + u ^ 2n + 2, and 3n + 3 < v < 4n + 2: 



A, 



:= -det 



C(l)| 



|C(3n + 2) 



We will deduce from the following lemma that D u<r ^ v is equal to zero unless u = j. 
2n + j + 2 and t> = 3n + j + 2, for 1 < j < n. 



(4.10) 
n + j + 1, or 



Lemma 4.5. Fix j € {1, 2, . . . , n} and call T,' n , the submatrix of obtained by deleting any two 
columns indexed by two elements of Ij . It holds that any 3n x 3n-minor of Y,' n is equal to zero. 



Proof. We will keep the notation R(i) from the proof of Proposition 14.41 We now prove the lemma 
first for j = 1, then j = n, and then finally for 1 < j < n. 

For the case j = 1, we focus on the reactions 1 <^ ra + 2 — >■ 2, 2n + 3 <^ 3ra + 3 — >■ 2n + 2, and 
3n + 4 -> 2n + 3. If we delete C(l) and C(n + 2), or C(2n + 3) and C(3n + 3), then the rows of 
corresponding to R(n + 2) or R(3n + 3) will be equal to zero and the minor will be zero. 

If we delete C(l) and C(2n+3) (or C(3n+3)), or we delete C{n + 2) and C(2n+3) (or C(3n+3)), 
the rows corresponding to R{n + 2) and R(3n + 3) will have only one entry different from zero and 
the determinant will be obviously zero if the column corresponding to any of this entries is not 
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considered, or it will be the product of two constants and a (3n — 2) x (3n — 2)-minor that does 
not include the columns C(l), C(n + 2), C(2n + 3), C(3n + 3) nor the rows R(n + 2), R(3n + 3). 

It is important to notice that the columns of A* carry the information of the reactions whose 
source (educt) is the corresponding complex, therefore, C{t) carries the information of the reaction 
whose source is the i-th. complex. As the only complexes that generate reactions whose product 
is the (n + 2)-nd or (3n + 3)-rd complexes are the first and (2n + 2) complexes, respectively, it 
follows that the columns that are being considered in this new (3n — 2) x (3n — 2)-minor carry 
the information of reactions that do not end in either the (n + 2)-nd or the (3n + 3)-rd complexes. 
Hence the sum of the rows in this new submatrix, and therefore the minor as well, is equal to zero. 

For j = n, the analysis is similar. 

For 1 < j < n we focus on the reactions j T± n+j+1 — > j+1 and 2n+j+2 ^ 3n+j+2 — > 2n+j+l. 
If we delete C(j) and C(n + j + 1), or C(2n + j + 2) and C(3n + j + 2), then the rows of T/ n 
corresponding to R(n + j + 1) or i?(3n + j + 2) will be equal to zero and the minor will be zero. 

If we delete C(j) and C(2n + j + 2) (or C(3n+j + 2)), or we delete C(n+j + 1) and C(2n+j + 2) 
(or C(3n+ j + 2)), the rows corresponding to R(n + j + 1) and R(3n+ j + 2) will have only one entry 
different from zero, and thus the determinant will be obviously zero if the column corresponding 
to any of these entries is not considered. Otherwise it will be the product of two nonzero rate 
constants and a (3ra — 2) x (3n — 2)-minor that does not include any of C(j), C(n + j + 1), C(2n + 
j + 2), C(3n + j + 2) nor any of R(n + j + l),R(3n +j + 2). 

But deleting these columns is equivalent to not considering the reactions whose sources (educts) 
are the complexes j,n + j + l,2n + j + 2, or 3n + j + 2. This disconnects the graph into four linkage 
classes, so this new graph gives a Laplacian matrix formed by four blocks. The rows of E n that we 
are considering in T,' n come from adding rows of the first and third blocks of A l K , or the second and 
fourth ones; and the last rows of E n , which correspond to intermediary species, clearly belong to 
only one of the blocks. Then, this new submatrix of T,' n can be reordered into a two-block matrix, 
for which the sums of the rows in each block are zero. Hence, the matrix obtained from Ti' n without 
these four columns and two rows has rank at most 3ra — 3 and therefore any (3n — 2) x (3n — 2)-minor 
will be zero. □ 

We are now ready to prove Theorem 14.31 



Proof of Theorem \4-3\ Due to Lemma 14.51 f° r a 3n x 3n-minor of T,' n to be different from zero, we 
must obtain these 3n columns by choosing three from each group indexed by Ij, for 1 < j < n. In 
fact, any 3n x 3n- minor of T,' n that includes three columns from each group of four indexed by Ij, 
for 1 < j < n, is always nonzero due to Proposition 14.41 

We can now solve system (|4.3|) by applying Cramer's rule. Recall the notation from (|4. 10[) : 



Vl 


-1 






~ IT 


_ v 3n _ 





<H-3n+3 



Di 



D 



3n+2«3n+3 



D 



3n+2-H>4n+2 



V3n+1 



V4n 



By Lemma 14.51 we already know that in the 3n x n-matrix in the right-hand side above, the only 
nonzero entries are Dj, D n+ j + i, and Z?2n+j+2- This gives us a description of ker(S n ), which has a 
basis of the following form: 

{e n+ i,e 2 „+2}U{6 1 ,6 2 ,...,6 n } 

for V as in (jlTTD . 

This proves that the n-site phosphorylation system satisfies Condition 13.11 for the partition 
Ji, I2, ■ ■ ■ , I-n+2 and the basis of ker(S n ), {b , b 2 , . . . ,b n , e n+ \, e2n+2}, described above. 

We now prove that the n-site phosphorylation system additionally satisfies Conditions l3.4l and l3.61 
Condition 13.41 is satisfied immediately by Proposition 14.41 With respect to Condition 13.61 we notice 
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that the subspace spanned by the columns of the matrix A has the following basis: 

{e2n+j+l — Cj — ^2n+j+l ~ Cn+j+1, €-2n+j+l ~ ~ £3n+3 | 1 < j < n } ■ (4-H) 

0. Hence, equation (13. 6p is trivially 



Therefore, the dimension of the image of A is 3n, so ker(A) 
satisfied, as noted in Remark 13.71 

Then, by Theorem 13.81 it is immediate that the n-site phosphorylation system has toric steady 
states that are positive and real. Finally, for a parametrization of the steady state locus, let us 
consider the following matrix: 



A 





1 



1 2 



1 1 



1 
1 1 




D 3x(3n+3) 



It has maximal rank, and its kernel equals the span of all the differences y,- 2 — j/,- a , for ji , j2 E Ij , 
where 1 < j < n + 2, shown in (|4.11j) . After applying Theorem 13.111 we are left to see that the 
point x defined in the statement of the present theorem is a positive steady state of the system. 
But it is easy to check that x is a positive steady state by applying Theorem 13.31 to the following 
binomials: 



DxjX 3n+2 - DjX2n 



+3+1 



Dx 



n+ j- 



-1 ~ Dn+j+lX2n+j+l, Dx j+1 X 3n+3 - D 2n +j+2X2n+j+l, for 1 < j < U. 



This completes the proof. 
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5. MULTISTATIONARITY FOR SYSTEMS WITH TORIC STEADY STATES 

In this section we focus on the capacity of a chemical reaction system with toric steady states 
to exhibit multiple steady states. Following prior work of Conradi et al. [5] and Holstein [19], we 
make use of an alternative notation for reaction systems to obtain a characterization of steady states 
(Proposition I5.2| ). This result is used to prove a criterion for the existence of multistationarity for 
systems with toric steady states that satisfy Conditions 13.11 13.41 and 13.61 (Theorem I5.5|) . At the 
end of this section, we make the connection to a related criterion of Feinberg. 

Often a chemical reaction system has a continuum of steady states, as long as one steady state 
exists. However, as defined earlier (and as it is in Chemical Engineering), multistationarity refers to 
the existence of multiple steady states within one and the same stoichiometric compatibility class. 
In general one is interested in situations where the steady state locus intersects a stoichiometric 
compatibility class in a finite number of points [13]. In Computational Biology one is sometimes 
interested in situations where the steady state locus intersects an affine subspace distinct from 
translates of the stoichiometric subspace S [15]. Here we define multistationarity with respect to 
a linear subspace in the following way. Consider a matrix Z G R SX|J , where q is a positive integer. 
We say that the chemical reaction system x = E ■ ^f(x) exhibits multistationarity with respect to the 
linear subspace ker(Z t ) if and only if there exist at least two distinct positive steady state vectors 
x , x 2 G suc h that their difference lies in ker(Z*); in other words the following equations must 
hold: 

^■^(x 1 ) = (5.1a) 
S-^(x 2 ) = (5.1b) 
Z l x x = Z f x 2 . (5.1c) 

Note that if the columns of Z form a basis for iS -1- , one recovers the usual definition of multista- 
tionarity given in Section 12.21 In this case, Equation (|5.1cj) states that the steady states x 1 and 
x 2 belong to the same stoichiometric compatibility class, and we simply speak of multistationarity, 
omitting the linear subspace we are referring to. 
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5.1. Second representation of a chemical reaction system. We now introduce a second 
representation of the differential equations that govern a chemical reaction system (|2.3j) ; this will 
prove useful for the characterization of steady states (Proposition I5.2p and for establishing the 
capacity of a chemical reaction network for multistationarity. Letting r denote the number of 
reactions of a chemical reaction network G, we fix an ordering of these r reactions and define the 
incidence matrix X G { — 1, 0, l} mxr of the network to be the matrix whose i-th column has a 1 in 
the row corresponding to the product complex of the i-th reaction and a —1 for the educt (reactant) 
complex. Then the (s x r)-matrix product 

N := Y l X (5.2) 

is known as the stoichiometric matrix. Thus, the i-th column of N is the reaction vector corre- 
sponding to reaction i. Next we define the educt- complex matrix 

y ■= [yi, 2/2, • • • , y r ] , (5-3) 

where the column yi of y is defined as the vector of the educt complex of the i-th reaction. Now 
we can define the vector of educt complex monomials 

4>{x) := X y\ x^Y . (5.4) 

We also define k G M> to be the vector of reaction rate constants: ki is the rate constant of the 
i-th reaction (that is, ki = Ki>ji where the i-th reaction is from the complex x Vi ' to x v i'). We now 
give a second formulation for a chemical reaction system (|2.3|) (cf. |16j): 

x = N diag(fc) <f>(x) . (5-5) 

Both formulations of a chemical reaction system given in equations (|2.3p and (|5.5p lead to the same 
system of ODEs and hence are equivalent. This can be made explicit by way of the doubling matrix 
D of dimension m x r which relates y and Y via y = Y D. Here the i-th column vector of D is 
defined as the unit vector e,- of M m such that yj is the educt (reactant) complex vector of the i-th 
reaction. From 

x = N diag(fc) <j){x) = Y l X diag(fc) D l (x) = m(x) , 
it follows that 4>{x) = D t ^f(x) and A l K = X diag(fc) D l . 

Example 5.1. For the 1-site phosphorylation network \2.J$ , one obtains the matrices 

1 " 



110 
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1 1 

1 " 

1 

110 

1 1 ' 

1 
1 

and the monomial vector <f>{x) = (x±X5, x$, xs, X2Xq, X4, x^f. 

It follows from the differential equations (|5.5|) that a positive concentration vector x £ R^g is a 
steady state for the chemical reaction system defined by the positive reaction rate constant vector 
k if and only if 

diag(A;) <t>{x) G ker(iV) nM^ . 
We now recognize that the set ker(iV) n M>Q' ^ nonempty, is the relative interior of the pointed 
polyhedral cone ker(iV) n!R>g. To utilize this cone, we collect a finite set of generators (also called 
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"extreme rays") of the cone ker(iV) n M> as columns of a non- negative matrix M. Up to scalar 
multiplication, generators of a cone are unique and form a finite set; as the cone of interest arises 
as the intersection of an orthant with a linear subspace, the generators are the vectors of the cone 
with minimal support with respect to inclusion. (Background on polyhedral cones can be found in 
the textbook of Rockafellar [27].) Letting p denote the number of generators of the cone, we can 
use M to express the condition for a positive vector x G M>q to be a steady state of the chemical 
reaction system in the following way: 

diag(fc) <f>(x) = MX , for some A G M| with MAG M. r >0 . (5.6) 

Note that this proves the following result which appears in [5]: 

Proposition 5.2 (Characterization of steady states of chemical reaction systems). For a chemical 
reaction network G, let M denote a corresponding generator matrix as defined above. Then a 
positive vector x G M> * s a steady state for the chemical reaction system defined by reaction rate 
vector k G M> , if and only if there exists a vector A G M> such that 

k = diag ((0(x)) _1 MX and MA€l r >0 . (5.7) 

We now note that outside of a degenerate case, any positive concentration vector can be a steady 
state for appropriately chosen rate constants k. 

Remark 5.3. We now comment on the degenerate case of a network for which the set ker(iV) nM> 
is empty. First, this case is equivalent to either of the following three conditions: (i) there is no 
positive dependence among the reaction vectors (yj — yi), (ii) the cone ker(iV) n M> is contained 
in a coordinate hyperplane, and (iii) the generator matrix M has at least one zero row. Now, 
in this degenerate case, it is clear that for any choice of reaction rate constants, the chemical 
reaction system has no positive steady states. This is because if x* G M> ^ s a steady state for 
the system with reaction rate constants n%j, then the numbers Oy := Kij ■ {x*) Vi witness to the 
positive dependence among the reaction vectors (yj — yi)'s. Outside of this degenerate case, it 
follows from Proposition 15.21 that there exists a vector of reaction rate constants k for which the 
resulting chemical reaction system has a positive steady state. Moreover, in this case any positive 
vector x can be a steady state, by choosing k as in equation (|5.7jl for some valid choice of A G K> - 

Using our new notation, we return to the question of existence of steady states. 

Remark 5.4. Recall the content of Corollary l3.12l for a chemical reaction network for which a single 
partition works to satisfy Condition 13.11 for all choices of positive rate constants, the set of rate 
constant vectors k that yield systems with positive steady states is the semialgebraic set of M> 
defined by Conditions 13.41 and 13.61 We now note that Proposition 15.21 implies that this set of rate 
constant vectors is the image of the following polynomial map: 

P : K s >0 xr ^ R r >Q 

(x,X) ^ diag(0(x)) _1 MX , 

where T := {AG K> [ MA G K>ol- I n case that Condition 13.11 holds and Condition 13.61 is 
trivially satisfied (i.e. A has full row rank), the image of (3 is cut out by the inequalities defined 
by Condition 13.41 

5.2. Main result on multistationarity. We now make use of Proposition 15.21 to examine which 
chemical reaction systems with toric steady states exhibit multistationarity. We first note that in 
the setting of Section O the set of differences lnx 1 — lnx 2 , where x 1 and x 2 are positive steady 
states for the same system, form a linear subspace. As before, the notation "In a;" for a vector 
x G R> denotes the vector (lnxi,lnx2, . . . ,lnx s ) G R s ; similarly we will make use of the notation 
"e x " to denote component- wise exponentiation. 
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Our next theorem, the main result of this section, is a consequence of Lemma 1]. It states 
that a network that satisfies Condition 13.11 has the capacity for multistationarity if and only if two 
subspaces, namely im(A') and S, both intersect non-trivially some (possibly lower-dimensional) 
orthant {x G M s | sign(x) = oj} defined by a sign vector u) G { — ,0, +} s . We remark that this is a 
matroidal condition. Related ideas appear in work of Feinberg |14j . and details on the connection 
between our work and Feinberg's appears at the end of this section. 

Theorem 5.5 (Multistationarity for networks with toric steady states). Fix a chemical reaction 
network G with s species and m complexes, and let Z G Z sx<? be an integer matrix, for some positive 
integer q. Assume that the cone ker(iV)nR>Q is not contained in any coordinate hyperplane. Assume 
moreover that there exists a partition I\, 1%, ■ ■ ■ , Id °f the m complexes of G such that Condition \3.1\ 
is satisfied for all rate constants. 

Recall the matrix A for this partition from the proof of Theorem \3.11\ Then there exists a reaction 
rate constant vector such that the resulting chemical reaction system exhibits multistationarity with 
respect to the linear subspace ker(Z*) if and only if there exists an orthant ofM s that both subspaces 
im(A') and ker (Z*) intersect nontrivially. More precisely, given nonzero vectors a G im(^4*) and 
a G ker (Z*) with 

sign(a) = sign(cr) , (5-8) 

then two steady states x 1 and x 2 and a reaction rate constant vector k that witness multistationarity 
(that is, that satisfy equations (15. lap . (|5.1bj) . and (15. left ) arise in the following way: 



(x\). 1 = Jsfe'tf*^ (5.9) 



where Xi denotes an arbitrary positive number, and 



.,2 



x 



diag(e°)x 1 (5.10) 



k = diag(</.(rE i ))" 1 MA , (5.11) 

for any non-negative vector A G M> for which MX G M^ . Conversely, any witness to mul- 
tistationarity with respect to ker (Z*) (given by some x 1 , x 2 G K^q, and k G ^-yo) arises from 
equations ()5.8[) . (|5.9[) . (|5.10|) . and (|5.1ip for some vectors a G im(yl*) and a G ker (Z*) that have 
the same sign. 

Proof. Assume that there exist nonzero vectors a G im(j4 4 ) and a G ker (Z*) having the same 
sign. First note that the vectors x , x 2 , and k defined by (|5.9|) . (|5.1U|) . and (|5.1ip . respectively, are 
positive because a and a have the same sign and because the cone ker(iV) n M> is not contained 
in a coordinate hyperplane. By Proposition I5.2| equation (|5,lip implies that x 1 is a steady state 
of the system defined by k. We now claim that x 2 too is a steady state of the same system. This 
follows from Theorem 13.111 because the difference between lnx 1 and Inx 2 is in im(^4'): 

lna^-lnz 2 = —a G im(A*) . 

Conversely, assume that vectors x , x 2 , and k are a witness to multistationarity with respect 
to ker(Z*). Let us now construct appropriate vectors a and a. By Theorem 13.111 the vector 
a := lnx 2 — lnx 1 is in im(^4*). Next, we define a G M s by o~i = (e ai — l)xj if «j 7^ and cxj = 
if on = 0, so by construction, a and a have the same sign. In addition, equations (|5.9p and (|5.10p 
easily follow for these values of a and a. We also see that 

—a = x 1 — x 2 G ker(Z*) , 

so a G ker(Z'). Finally, Proposition ^ . 2 l implies that there exists a valid A G M> that satisfies (|5.1ip . 

□ 
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Remark 5.6. If a chemical reaction system defined by reaction rate constants k* and a partition 
of its complexes satisfy Conditions 13.14 13.41 and 13.61 (but not necessarily for other choices of rate 
constants) , then the equations (|5.8p , (|5.9p , (|5.10p , and (|5.1ip in Theorem 15.51 still characterize 
multistationarity. In other words, x 1 and x 2 are two steady states that demonstrate that the 
system defined by k* has the capacity for multistationarity with respect to ker(Z*) if and only if 
there exist a £ im(A*), a G ker(Z*), and A S R> such that those four equations hold. 

Example 5.7 (Triangle network, continued). We return to the Triangle network analyzed in Ex- 
amples 0O1 and \3. 1 (A The stoichiometric subspace is 

ker(S) = S = span{(l,-l)} . 

In the toric setting (recall that this is when K31 = the partition for which the system satisfies 
Condition \3.1\ is {1, 2}, {3} ; so a matrix A for which 

kei(A) = span{y 2 -yi} = span{(2,-2)} 

is A = [1 1]. We can see that the subspaces ker(Z') and im(A*) = span{(l, 1)} do not both 
intersect some orthant nontrivially. So Theorem 15.51 allows us to conclude that no system (for 
which K31 = arising from the Triangle network exhibits multistationarity. 

Although the capacity of the Triangle network to exhibit multistationarity is easily determined 
directly, without the need to apply Theorem 15.51 it is more difficult in the case of the multisite 
phosphorylation system. Recall that we proved in Theorem 14.31 that any n-site phosphorylation 
system satisfies Condition 13.11 with the same partition (for fixed n). Hence, Theorem 15.51 can be 
used to compute the semialgebraic set of reaction rate constants k that give rise to multistationarity 
for the phosphorylation networks. This was performed by Conradi et al. (for the 2-site network) [5] 
and Holstein (for the general n-site network) [19]; multistationarity is possible only for n > 2. 
Results on the number of steady states of phosphorylation systems appeared in work of Wang and 
Sontag [36] and is the focus of a forthcoming work of the authors [1] . 

5.3. Connection to related results on multistationarity. We now make the connection be- 
tween our results on the capacity of a chemical reaction network to exhibit multistationarity and 
related results of Feinberg |14j . A regular network is a network for which (i) ker(iV) n M> 7^ 0, 
(ii) each linkage class contains a unique terminal strong linkage class, and (iii) removing the re- 
action^) between any two adjacent complexes in a terminal strong linkage class disconnects the 
corresponding linkage class. Recall from Remark 15.31 that condition (i) in this definition is simply 
the requirement that the reaction vectors yj — yi are positively dependent, and that this condition 
is necessary for the existence of positive steady states. Recall that the deficiency of a chemical 
reaction network was discussed in § 12.31 

We now can explain the relationship between Feinberg's result and ours. Feinberg examined 
regular deficiency-one networks, while we are concerned with networks for which there exists a 
partition that satisfies Condition 13.61 (for all rate constants). In these respective settings, both 
Theorem 4.1 and Corollary 4.1 of |14j and Theorem 15.51 in this article state that a certain subset 
of R s and the stoichiometric subspace both intersect the same orthant non-trivially if and only if 
the network has the capacity for multistationarity. In the result of Feinberg, this set is a union 
of certain polyhedral cones, while in our case, this set is the image of A 1 . In both cases, this set 
consists of all vectors ln(c*/c**), where c* and c** are steady states arising from the same rate 
constants. As an illustration, see Example 15.81 below. 

Let us now explain how the two results are complementary. First, there are some networks for 
which only Feinberg's results apply. For example, consider any network for which the union of 
polyhedral cones obtained from Feinberg's results is not a linear space. Additionally, for some 
networks, only our results apply. As an example, the n > 1 multisite networks have deficiency 



28 



MERCEDES PEREZ MILLAN, ALICIA DICKENSTEIN, ANNE SHIU, AND CARSTEN CONRADI 



greater than one. Finally, for some networks, both our results and Feinberg's apply, such as in the 
following example. 

Example 5.8. The 1-site phosphorylation network of Example \2.1\ is regular and has deficiency 
one. In this case, both the image of A 1 and Feinberg's union of cones are the subspace o/M 6 spanned 
by the three vectors (ei + ei + e% + e^), [e-i + e^ + e^ + e§), and (e^ + e^ + es + e%). So in this 
instance, our Theorem \5.5\ and Feinberg's Corollary 4-1 of [14] coincide. 

Finally we note that the proofs of both results make use of special structure of ker(S). In our 
case, we assume the existence of a basis with disjoint support. For Feinberg's results, there is 
a non-negative basis where the supports of the first L basis vectors correspond exactly to the L 
terminal strong linkage classes, and the last basis vector is the all-ones vector (here L denotes the 
number of terminal strong linkage classes) . 
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